library(tidyverse) ; library(reshape2) ; library(glue) ; library(plotly) ; library(dendextend)
library(RColorBrewer) ; library(viridis) ; require(gridExtra) ; library(GGally)
library(expss)
library(polycor)
library(foreach) ; library(doParallel)
suppressWarnings(suppressMessages(library(WGCNA)))
SFARI_colour_hue = function(r) {
pal = c('#FF7631','#FFB100','#E8E328','#8CC83F','#62CCA6','#59B9C9','#b3b3b3','#808080','gray','#d9d9d9')[r]
}
Load preprocessed dataset (preprocessing code in 19_11_14_data_preprocessing.Rmd) and clustering (pipeline in 20_01_23_WGCNA.Rmd)
# Gupta dataset
load('./../Data/preprocessed_data.RData')
datExpr = datExpr %>% data.frame
DE_info = DE_info %>% data.frame
# GO Neuronal annotations
GO_annotations = read.csv('./../../FirstYearReview/Data/GO_annotations/genes_GO_annotations.csv')
GO_neuronal = GO_annotations %>% filter(grepl('neuron', go_term)) %>%
mutate('ID'=as.character(ensembl_gene_id)) %>%
dplyr::select(-ensembl_gene_id) %>% distinct(ID) %>%
mutate('Neuronal'=1)
# SFARI Genes
SFARI_genes = read_csv('./../../../PhD-Models/FirstPUModel/Data/SFARI/SFARI_genes_08-29-2019_with_ensembl_IDs.csv')
SFARI_genes = SFARI_genes[!duplicated(SFARI_genes$ID) & !is.na(SFARI_genes$ID),]
# Clusterings
clusterings = read_csv('./../Data/clusters.csv')
# Update DE_info with SFARI and Neuronal information
genes_info = DE_info %>% mutate('ID'=rownames(.)) %>% left_join(SFARI_genes, by='ID') %>%
mutate(`gene-score`=ifelse(is.na(`gene-score`), 'None', `gene-score`)) %>%
left_join(GO_neuronal, by='ID') %>% left_join(clusterings, by='ID') %>%
mutate(Neuronal=ifelse(is.na(Neuronal), 0, Neuronal)) %>%
mutate(gene.score=ifelse(`gene-score`=='None' & Neuronal==1, 'Neuronal', `gene-score`),
significant=padj<0.05 & !is.na(padj))
rm(DE_info, GO_annotations, clusterings)
print(paste0('Dynamic Tree leaves ', sum(genes_info$DynamicTree=='gray'), ' genes without cluster (',
round(mean(genes_info$DynamicTree=='gray')*100), '%)'))
## [1] "Dynamic Tree leaves 12639 genes without cluster (84%)"
print(paste0('Dynamic Hybrid leaves ', sum(genes_info$DynamicHybrid=='gray'), ' genes without cluster (',
round(mean(genes_info$DynamicHybrid=='gray')*100), '%)'))
## [1] "Dynamic Hybrid leaves 520 genes without cluster (3%)"
Dynamic Tree leaves many more genes without a cluster, so perhaps this is not the best option with this new dataset…
There seems to be a relation between DE and module membership, being DE a more restrictive condition than being assigned to a cluster. In both cases, most of the genes were left out.
This could be a consequence of the sva batch correction, because all the patterns that survived were diagnosis related, so it could have left all the genes with no relation to ASD without any recognisable pattern, which left them at the top of the dendrogram, which are the ones the clustering algorithm didn’t assign to any module
pca = datExpr %>% prcomp
plot_data = data.frame('ID'=rownames(datExpr), 'PC1' = pca$x[,1], 'PC2' = pca$x[,2]) %>%
left_join(genes_info, by='ID') %>% mutate('hasCluster'=DynamicTree!='gray',
'hasSFARIScore'=`gene-score`!='None') %>%
apply_labels(`gene-score`='SFARI Gene score', DynamicTree = 'Dynamic Tree Algorithm',
significant = 'Differentially Expressed', hasCluster = 'Belongs to a Module',
hasSFARIScore = 'Has a SFARI Score', syndromic = 'Has syndromic tag')
p1 = plot_data %>% ggplot(aes(PC1, PC2, color=hasCluster)) + geom_point(alpha=0.2) +
theme_minimal() + ggtitle('Genes are assigned to a cluster') + theme(legend.position='bottom')
p2 = plot_data %>% ggplot(aes(PC1, PC2, color=significant)) + geom_point(alpha=0.2) +
theme_minimal() + ggtitle('Genes were found to be DE') + theme(legend.position='bottom')
grid.arrange(p1, p2, nrow=1)
rm(pca, p1, p2)
Most of the genes that don’t have a cluster (98%) are not differentially expressed.
cro(plot_data$significant, list(plot_data$hasCluster, total()))
|  Belongs to a Module |  #Total | |||
|---|---|---|---|---|
| Â FALSEÂ | Â TRUEÂ | Â | ||
|  Differentially Expressed | ||||
| Â Â Â FALSEÂ | 12339 | 2121 | Â | 14460 |
| Â Â Â TRUEÂ | 300 | 321 | Â | 621 |
|    #Total cases | 12639 | 2442 |  | 15081 |
But we lose most of the SFARI genes if we do that
cro(plot_data$hasSFARIScore, list(plot_data$hasCluster, total()))
|  Belongs to a Module |  #Total | |||
|---|---|---|---|---|
| Â FALSEÂ | Â TRUEÂ | Â | ||
|  Has a SFARI Score | ||||
| Â Â Â FALSEÂ | 11950 | 2269 | Â | 14219 |
| Â Â Â TRUEÂ | 689 | 173 | Â | 862 |
|    #Total cases | 12639 | 2442 |  | 15081 |
print(paste0(sum(plot_data$hasSFARIScore & !plot_data$hasCluster), ' of the SFARI genes (',
round(100*sum(plot_data$hasSFARIScore & !plot_data$hasCluster)/sum(plot_data$hasSFARIScore)),
'%) are not assigned to any cluster'))
## [1] "689 of the SFARI genes (80%) are not assigned to any cluster"
Separated by scores
cro(plot_data$`gene-score`, list(plot_data$hasCluster, total()))
|  Belongs to a Module |  #Total | |||
|---|---|---|---|---|
| Â FALSEÂ | Â TRUEÂ | Â | ||
|  SFARI Gene score | ||||
| Â Â Â 1Â | 15 | 9 | Â | 24 |
| Â Â Â 2Â | 53 | 11 | Â | 64 |
| Â Â Â 3Â | 153 | 29 | Â | 182 |
| Â Â Â 4Â | 335 | 79 | Â | 414 |
| Â Â Â 5Â | 114 | 43 | Â | 157 |
| Â Â Â 6Â | 19 | 2 | Â | 21 |
|    None | 11950 | 2269 |  | 14219 |
|    #Total cases | 12639 | 2442 |  | 15081 |
rm(plot_data)
The main difference between algorithms is that Dynamic Hybrid clusters outlier genes and Dynamic Tree leaves them out, so Dynamic Tree would give me a ‘cleaner’ group of genes to work with but losing too many genes, and Dynamic Hybrid would give me less and more balanced clusters, but the Dynamic Tree algorithm leaves too many genes (including most of the SFARI genes) without a cluster, so I’m going to use the Dynamic Hybrid results.
Using the clustering from Dynamic Hybrid
clustering_selected = 'DynamicHybrid'
genes_info$Module = genes_info[,clustering_selected]
*The colour of the modules is the arbitrary one assigned during the WGCNA algorithm, where the gray cluster actually represents all the genes that were left without a cluster (so it’s not actually a cluster).
cat(paste0('The Dynamic Hybrid algorithm created ', length(unique(genes_info$Module))-1, ' modules and leaves ',
sum(genes_info$Module=='gray'), ' genes without a module.\n'))
## The Dynamic Hybrid algorithm created 123 modules and leaves 520 genes without a module.
table(genes_info$Module)
##
## #00A4FF #00A7FF #00A9FF #00ABFC #00AEFA #00B0F7 #00B1F4 #00B3F1 #00B5ED
## 29 621 157 36 147 47 136 70 201
## #00B7EA #00B817 #00B8E6 #00B92B #00B9E2 #00BA38 #00BADE #00BB44 #00BB4D
## 44 61 28 136 270 312 65 303 18
## #00BC56 #00BCD9 #00BD5F #00BD66 #00BDD0 #00BDD5 #00BE6E #00BECC #00BF75
## 54 319 65 89 27 31 84 57 83
## #00BF7C #00BFC2 #00BFC7 #00C083 #00C089 #00C08F #00C096 #00C0B2 #00C0B8
## 308 33 104 363 67 127 83 51 129
## #00C0BD #00C19C #00C1A1 #00C1A7 #00C1AD #20B700 #37A1FF #38B600 #48B500
## 74 44 29 247 30 85 156 113 20
## #4F9FFF #54B400 #5FB200 #619CFF #69B100 #7099FF #72B000 #7AAF00 #7D96FF
## 49 177 61 180 111 88 33 140 114
## #81AD00 #88AC00 #8993FF #8FAA00 #9490FF #95A900 #9BA800 #9E8DFF #A1A600
## 677 60 61 51 78 25 130 90 190
## #A6A400 #A78AFF #ABA300 #AF86FF #B0A100 #B5A000 #B783FF #B99E00 #BE80FF
## 280 15 32 275 135 70 274 90 244
## #BE9C00 #C29A00 #C57DFF #C69900 #CA9700 #CC7AFF #CE9500 #D277FF #D29300
## 13 38 226 36 84 55 87 45 146
## #D59100 #D774FD #D98F00 #DC71FA #DC8D00 #DF8B00 #E16FF7 #E28900 #E58706
## 363 82 97 33 57 26 116 37 228
## #E66CF4 #E88523 #EA6AF1 #EB8333 #ED68ED #ED8140 #F066E9 #F07F4A #F27D54
## 216 38 35 187 302 87 44 77 80
## #F365E6 #F47A5D #F663E1 #F67865 #F862DD #F8766D #FA7474 #FB61D9 #FB727C
## 106 48 137 34 33 27 85 181 92
## #FC61D4 #FD7083 #FE61CF #FE6E8A #FF61C5 #FF61CA #FF62BA #FF62C0 #FF63B5
## 44 71 17 281 108 65 117 268 92
## #FF64AF #FF66A9 #FF67A3 #FF699D #FF6B97 #FF6C90 gray
## 53 138 149 83 183 31 520
plot_data = table(genes_info$Module) %>% data.frame %>% arrange(desc(Freq))
ggplotly(plot_data %>% ggplot(aes(x=reorder(Var1, -Freq), y=Freq)) + geom_bar(stat='identity', fill=plot_data$Var1) +
ggtitle('Module size') + ylab('Number of genes') + xlab('Module') + theme_minimal() +
theme(axis.text.x = element_text(angle = 90)))
In the WGCNA documentation they use Pearson correlation to calculate correlations, I think all of their variables were continuous. Since I have categorical variables I’m going to use the hetcor function, that calculates Pearson, polyserial or polychoric correlations depending on the type of variables involved.
I’m not sure how the corPvalueStudent function calculates the p-values and I cannot find any documentation…
Compared correlations using Pearson correlation and with hetcor and they are very similar, but a bit more extreme with hetcor. The same thing happens with the p-values.
datTraits = datMeta %>% dplyr::select(Diagnosis, brain_lobe, Gender, Age, PMI, SiteHM) %>% rename('Batch' = SiteHM)
# Recalculate MEs with color labels
ME_object = datExpr %>% t %>% moduleEigengenes(colors = genes_info$Module)
MEs = orderMEs(ME_object$eigengenes)
# Calculate correlation between eigengenes and the traits and their p-values
# I was originally using ML=TRUE in the hetcor function, but in the Gandal dataset that created an NA and in this dataset it took too long to run
moduleTraitCor = MEs %>% apply(2, function(x) hetcor(x, datTraits)$correlations[1,-1]) %>% t
rownames(moduleTraitCor) = colnames(MEs)
colnames(moduleTraitCor) = colnames(datTraits)
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nrow(datExpr))
# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)
# In case there are any NAs
if(sum(!complete.cases(moduleTraitCor))>0){
print(paste0(sum(is.na(moduleTraitCor)),' correlation(s) could not be calculated'))
}
rm(ME_object)
Modules have strong correlations with Diagnosis (although not as strong as with Gandal’s dataset) with really small p-values. Because of the batch effect, there are strong correlations to processing site as well, although they aren’t as strong as with Diagonsis, and there’s also one module with a strong correlation with gender.
It’s a good sign that the gray module has a low correlation with diagnosis, since we know its composed mainly of not differentially expressed genes.
# Sort moduleTraitCor by Diagnosis
moduleTraitCor = moduleTraitCor[order(moduleTraitCor[,1], decreasing=TRUE),]
moduleTraitPvalue = moduleTraitPvalue[order(moduleTraitCor[,1], decreasing=TRUE),]
# Create text matrix for the Heatmap
textMatrix = paste0(signif(moduleTraitCor, 2), ' (', signif(moduleTraitPvalue, 1), ')')
dim(textMatrix) = dim(moduleTraitCor)
labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits), yLabels = gsub('ME','',rownames(moduleTraitCor)),
yColorWidth=0, colors = brewer.pal(11,'PiYG'), bg.lab.y = gsub('ME','',rownames(moduleTraitCor)),
textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.8, cex.lab.y = 0.75, zlim = c(-1,1),
main = paste('Module-Trait relationships'))
diagnosis_cor = data.frame('Module' = gsub('ME','',rownames(moduleTraitCor)),
'MTcor' = moduleTraitCor[,'Diagnosis'],
'MTpval' = moduleTraitPvalue[,'Diagnosis'])
genes_info = genes_info %>% left_join(diagnosis_cor, by='Module')
## Warning: Column `Module` joining character vector and factor, coercing into
## character vector
rm(moduleTraitCor, moduleTraitPvalue, datTraits, textMatrix, diagnosis_cor)
Modules with a high Module-Diagnosis correlation should have a high content of differentially expressed genes.
This is specially strong for modules with a negative correlation. Many modules have no differentially expressed genes because in total there aren’t that many of these genes, usually they are small modules (which can be seen by the size of the point in the plot).
plot_data = genes_info %>% group_by(Module, MTcor) %>% summarise(p = 100*mean(significant), size=n())
ggplotly(plot_data %>% ggplot(aes(MTcor, p)) + geom_hline(yintercept=mean(plot_data$p), color='gray', linetype='dotted') +
geom_point(color=plot_data$Module, alpha=0.8, aes(id=Module, size=size)) + theme_minimal() +
xlab('Modules ordered by Module-Diagnosis correlation') + ylab('Percentage of differentially expressed genes'))
Gene significance: is the value between the correlation between the gene and the trait we are interested in. A positive gene significance means the gene is overexpressed and a negative value means its underexpressed. (The term ‘significance’ is not very acurate because it’s not actually measuring statistical significance, it’s just a correlation, but that’s how they call it in WGCNA…)
Module Membership is the correlation of the module’s eigengene and the expression profile of a gene. The higher the Module Membership, the more similar the gene is to the genes that constitute the module. (I won’t use this metric yet)
# It's more efficient to iterate the correlations one by one, otherwise it calculates correlations between the eigengenes and also between the genes, which we don't need
# Check if MM information already exists and if not, calculate it
if(file.exists(paste0('./../Data/dataset_', clustering_selected, '.csv'))){
dataset = read.csv(paste0('./../Data/dataset_', clustering_selected, '.csv'))
dataset$Module = dataset[,clustering_selected]
} else {
############# 1. Calculate Gene Significance
GS_info = data.frame('ID' = rownames(datExpr),
'GS' = datExpr %>% apply(1, function(x) hetcor(x, datMeta$Diagnosis)$correlations[1,2])) %>%
mutate('GSpval' = corPvalueStudent(GS, ncol(datExpr)))
############# 2. Calculate Module Membership
#setup parallel backend to use many processors
cores = detectCores()
cl = makeCluster(cores-1)
registerDoParallel(cl)
# Create matrix with MM by gene
MM = foreach(i=1:nrow(datExpr), .combine=rbind) %dopar% {
library(polycor)
tempMatrix = apply(MEs, 2, function(x) hetcor(as.numeric(datExpr[i,]), x)$correlations[1,2])
tempMatrix
}
# Stop clusters
stopCluster(cl)
rownames(MM) = rownames(datExpr)
colnames(MM) = paste0('MM',gsub('ME','',colnames(MEs)))
# Calculate p-values
MMpval = MM %>% corPvalueStudent(ncol(datExpr)) %>% as.data.frame
colnames(MMpval) = paste0('MMpval', gsub('ME','',colnames(MEs)))
MM = MM %>% as.data.frame %>% mutate(ID = rownames(.))
MMpval = MMpval %>% as.data.frame %>% mutate(ID = rownames(.))
# Join and save results
dataset = genes_info %>% dplyr::select(ID, `gene-score`, clustering_selected, MTcor, MTpval) %>%
left_join(GS_info, by='ID') %>%
left_join(MM, by='ID') %>%
left_join(MMpval, by='ID')
write.csv(dataset, file = paste0('./../Data/dataset_', clustering_selected, '.csv'), row.names = FALSE)
rm(cores, cl)
}
Gene significance and Log Fold Chance are two different ways to measure the same thing, so there should be a concordance between them
But both variables agree with each other quite well
plot_data = dataset %>% dplyr::select(ID, MTcor, GS) %>% left_join(genes_info %>% dplyr::select(ID, gene.score), by='ID') %>%
left_join(genes_info %>% dplyr::select(ID, baseMean, log2FoldChange, significant, Module), by='ID') %>%
left_join(data.frame(MTcor=unique(dataset$MTcor)) %>% arrange(by=MTcor) %>%
mutate(order=1:length(unique(dataset$MTcor))), by='MTcor')
ggplotly(plot_data %>% ggplot(aes(GS, log2FoldChange)) + geom_point(color=plot_data$Module, alpha=0.2) +
geom_smooth(color='gray', se=FALSE) + theme_minimal() + xlab('Gene Significance') +
ggtitle(paste0('Correlation = ', round(cor(plot_data$log2FoldChange, plot_data$GS)[1], 2))))
In general, modules with the highest Module-Diagnosis correlation should have genes with high Gene Significance
Note: For the Module-Diagnosis plots, if you do boxplots, you lose the exact module-diagnosis correlation and you only keep the order, so I decided to compensate this downside with a second plot, where each point is plotted individually using their module’s Module-Diagnosis correlation as the x axis. I think the boxplot plot is easier to understand but the second plot contains more information, so I don’t know which one is better.
ggplotly(plot_data %>% ggplot(aes(order, GS, group=order)) + geom_hline(yintercept=0, color='gray', linetype='dotted') +
geom_boxplot(fill=unique(plot_data$Module)) + theme_minimal() +
xlab('Modules ordered by Module-Diagnosis correlation') + ylab('Gene Significance'))
ggplotly(plot_data %>% ggplot(aes(MTcor, GS)) + geom_hline(yintercept=0, color='gray', linetype='dotted') +
geom_point(color=plot_data$Module, alpha=0.1, aes(id=ID)) + geom_smooth(color='gray', alpha=0.3) +
theme_minimal() + xlab('Module-Diagnosis correlation') + ylab('Gene Significance'))
The same should happen with the Log Fold Change
ggplotly(plot_data %>% ggplot(aes(order, log2FoldChange, group=order)) + geom_hline(yintercept=0, color='gray', linetype='dotted') +
geom_boxplot(fill=unique(plot_data$Module)) +
theme_minimal() + xlab('Modules ordered by Module-Diagnosis correlation') + ylab('log2FoldChange'))
ggplotly(plot_data %>% ggplot(aes(MTcor, log2FoldChange)) + geom_hline(yintercept=0, color='gray', linetype='dotted') +
geom_point(color=plot_data$Module, alpha=0.1, aes(id=ID)) + geom_smooth(color='gray', alpha=0.3) +
theme_minimal() + xlab('Module-Diagnosis correlation') + ylab('log2FoldChange'))
In theory, there shouldn’t be a relation between module-diagnosis and mean expression, but in the the exploratory analysis, we saw that the overexpressed genes tended to have lower levels of expression than the overexpressed genes, and this pattern can be seen in these plots where the modules with negative Module-Diagonsis correlation have slightly higher levels of expression than the modules with positive Module-Diagnosis correlation, although this pattern is note very strong and all modules have similar levels of expression.
ggplotly(plot_data %>% ggplot(aes(order, log2(baseMean+1), group=order)) +
geom_hline(yintercept=mean(log2(plot_data$baseMean+1)), color='gray', linetype='dotted') +
geom_boxplot(fill=unique(plot_data$Module)) + theme_minimal() +
xlab('Modules ordered by Module-Diagnosis correlation') + ylab('log2(Mean Expression+1)'))
ggplotly(plot_data %>% ggplot(aes(MTcor, log2(baseMean+1))) + geom_point(alpha=0.2, color=plot_data$Module, aes(id=ID)) +
geom_hline(yintercept=mean(log2(plot_data$baseMean+1)), color='gray', linetype='dotted') + ylab('log2(Mean Expression+1)') +
geom_smooth(color='gray', alpha=0.3) + theme_minimal() + xlab('Module-Diagnosis correlation'))
All of the variables seem to agree with each other, Modules with a high correlation with Diagnosis tend to have genes with high values of Log Fold Change as well as high values of Gene Significance, and the gray module, which groups all the genes that weren’t assigned to any cluster tends to have a very poor performance in all of the metrics.
Since SFARI scores genes depending on the strength of the evidence linking it to the development of autism, in theory, there should be some concordance between the metrics we have been studying above and these scores…
SFARI scores 3 to 5 have a lower median than all genes that have a neuronal-related annotation <span style="color:‘red’>!
The group with the highest Gene Significance is SFARI score 6, which is supposed to be the one with the least amount of evidence suggesting a relation to autism <span style="color:‘red’>!
SFARI score 5 is the group with the lowest Gene Significance, with the same median than the genes without any type of Neuronal annotation <span style="color:‘red’>!
Neuronal annotated genes have higher Gene Significance than genes without any neuronal-related annotation (makes sense)
ggplotly(plot_data %>% ggplot(aes(gene.score, abs(GS), fill=gene.score)) + geom_boxplot() +
scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
ylab('abs(Gene Significance)') + xlab('SFARI Scores') + theme(legend.position='none'))
ggplotly(plot_data %>% ggplot(aes(gene.score, abs(MTcor), fill=gene.score)) + geom_boxplot() +
scale_fill_manual(values=SFARI_colour_hue(r=c(1:6,8,7))) + theme_minimal() +
ylab('abs(Module-Trait Correlation)') + xlab('SFARI Scores') + theme(legend.position='none'))
Not only are SFARI genes not consistent with the other measurements, but they seem to contradict them, although these differences are not as strong as with Gandal’s dataset.
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-redhat-linux-gnu (64-bit)
## Running under: Scientific Linux 7.6 (Nitrogen)
##
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
##
## locale:
## [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
## [5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
## [7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] WGCNA_1.68 fastcluster_1.1.25 dynamicTreeCut_1.63-1
## [4] doParallel_1.0.15 iterators_1.0.12 foreach_1.4.7
## [7] polycor_0.7-10 expss_0.10.1 GGally_1.4.0
## [10] gridExtra_2.3 viridis_0.5.1 viridisLite_0.3.0
## [13] RColorBrewer_1.1-2 dendextend_1.13.2 plotly_4.9.1
## [16] glue_1.3.1 reshape2_1.4.3 forcats_0.4.0
## [19] stringr_1.4.0 dplyr_0.8.3 purrr_0.3.3
## [22] readr_1.3.1 tidyr_1.0.0 tibble_2.1.3
## [25] ggplot2_3.2.1 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.1.5
## [3] Hmisc_4.2-0 plyr_1.8.5
## [5] lazyeval_0.2.2 splines_3.6.0
## [7] crosstalk_1.0.0 BiocParallel_1.20.1
## [9] GenomeInfoDb_1.22.0 robust_0.4-18.2
## [11] digest_0.6.23 htmltools_0.4.0
## [13] GO.db_3.10.0 fansi_0.4.1
## [15] magrittr_1.5 checkmate_1.9.4
## [17] memoise_1.1.0 fit.models_0.5-14
## [19] cluster_2.0.8 annotate_1.64.0
## [21] modelr_0.1.5 matrixStats_0.55.0
## [23] colorspace_1.4-1 blob_1.2.0
## [25] rvest_0.3.5 rrcov_1.4-7
## [27] haven_2.2.0 xfun_0.8
## [29] crayon_1.3.4 RCurl_1.95-4.12
## [31] jsonlite_1.6 genefilter_1.68.0
## [33] zeallot_0.1.0 impute_1.60.0
## [35] survival_2.44-1.1 gtable_0.3.0
## [37] zlibbioc_1.32.0 XVector_0.26.0
## [39] DelayedArray_0.12.2 BiocGenerics_0.32.0
## [41] DEoptimR_1.0-8 scales_1.1.0
## [43] mvtnorm_1.0-11 DBI_1.1.0
## [45] Rcpp_1.0.3 xtable_1.8-4
## [47] htmlTable_1.13.1 foreign_0.8-71
## [49] bit_1.1-15.1 preprocessCore_1.48.0
## [51] Formula_1.2-3 stats4_3.6.0
## [53] htmlwidgets_1.5.1 httr_1.4.1
## [55] ellipsis_0.3.0 acepack_1.4.1
## [57] pkgconfig_2.0.3 reshape_0.8.8
## [59] XML_3.98-1.20 farver_2.0.3
## [61] nnet_7.3-12 dbplyr_1.4.2
## [63] locfit_1.5-9.1 later_1.0.0
## [65] tidyselect_0.2.5 labeling_0.3
## [67] rlang_0.4.2 AnnotationDbi_1.48.0
## [69] munsell_0.5.0 cellranger_1.1.0
## [71] tools_3.6.0 cli_2.0.1
## [73] generics_0.0.2 RSQLite_2.2.0
## [75] broom_0.5.3 fastmap_1.0.1
## [77] evaluate_0.14 yaml_2.2.0
## [79] knitr_1.24 bit64_0.9-7
## [81] fs_1.3.1 robustbase_0.93-5
## [83] nlme_3.1-139 mime_0.8
## [85] xml2_1.2.2 compiler_3.6.0
## [87] rstudioapi_0.10 reprex_0.3.0
## [89] geneplotter_1.64.0 pcaPP_1.9-73
## [91] stringi_1.4.5 lattice_0.20-38
## [93] Matrix_1.2-17 vctrs_0.2.1
## [95] pillar_1.4.3 lifecycle_0.1.0
## [97] data.table_1.12.8 bitops_1.0-6
## [99] httpuv_1.5.2 GenomicRanges_1.38.0
## [101] R6_2.4.1 latticeExtra_0.6-28
## [103] promises_1.1.0 IRanges_2.20.2
## [105] codetools_0.2-16 MASS_7.3-51.4
## [107] assertthat_0.2.1 SummarizedExperiment_1.16.1
## [109] DESeq2_1.26.0 withr_2.1.2
## [111] S4Vectors_0.24.2 GenomeInfoDbData_1.2.2
## [113] mgcv_1.8-28 hms_0.5.3
## [115] grid_3.6.0 rpart_4.1-15
## [117] rmarkdown_1.14 Cairo_1.5-10
## [119] shiny_1.4.0 Biobase_2.46.0
## [121] lubridate_1.7.4 base64enc_0.1-3